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Abstract 

We investigate cascades of isochronous pitchfork bifurcations of straight-hne hbrating orbits in some 
two-dimensional Hamiltonian systems with mixed phase space. We show that the new bifurcated 
orbits, which are responsible for the onset of chaos, are given analytically by the periodic solutions 
of the Lame equation as classified in 1940 by Ince. In Hamiltonians with C2v symmetry, they 
occur alternatingly as Lame functions of period 2K and 4K, respectively, where 4K is the period 
of the Jacobi elliptic function appearing in the Lame equation. We also show that the two pairs of 
orbits created at period-doubling bifurcations of island-chain type are given by two different linear 
combinations of algebraic Lame functions with period 8K. 



1 Introduction 

One of the well-established routes to chaos in maps is the so-called Feigenbaum scenario which 
consists in a cascade of successive period-doubling bifurcations of pitchfork type. They were first 
discussed for the 1-dimensional logistic map by Feigenbaum [|[| and then also found in the area- 
conserving two-dimensional Henon map |^, |3| , although the numerical scaling constants found there 
differ from those in the 1-dimensional case. One of us (M.B.) has recently investigated Q similar 
cascades of pitchfork bifurcations occurring in 2-dimensional Hamiltonian systems with mixed dy- 
namics, whereby the scaling constants can be determined analytically and depend on the potential 
parameters. In the present paper, we shall show that the new orbits born at its bifurcations are, 
near the bifurcations points, analytically given by periodic solutions of a linear second-order differ- 
ential equation studied over 160 years ago by G. Lame |p, and therefore called the "periodic Lame 
functions" |^. They have been classified uniquely in 1940 by Ince who also derived their Fourier 
series expansions Q. We find that these not only reproduce accurately the periodic orbits found 
numerically by solving the equations of motion at the bifurcations, but in the Henon-Heiles and 
similar potentials the Lame functions can also be used to describe the evolution of the bifurcated 
orbits at higher energies. A particularly interesting case is the homogeneous quartic oscillator for 
which the Lame functions become finite polynomials in terms of Jacobi elliptic functions. Here 
we can also find analytical expressions for the algebraic Lame functions which describe the orbits 
created at period-doubling bifurcations of island-chain type. 



1 



2 Bifurcations of a straight-line librating orbit 



We start from an autonomous two-dimensional Hamiltonian of a particle with unit mass in a smooth 
potential V{x, y) 

H=\(pl+p'^+V{x,y). (1) 

Assume that there exists a straight-line librating orbit, called A, along the y axis, so that 

XA{t) = G, yA{t) = yA{t + TA) (2) 

are solutions of the equations of motion and Ta is the period of the A orbit. Its stability is 
obtained from the stability matrix Ma that describes the propagation of the linearized flow of a 
small perturbation 5x{t), Spx{t) = 5x[t) transverse to the orbit A: 

[ Sp.. iW J = (0) J ■ 

When —2 < trMA < +2, the orbit is stable, for |trMA| > 2 it is unstable. Marginally stable orbits 
with tr Ma = +2 occur in systems with continuous symmetries; in two dimensions this would imply 
integrability. We investigate here only non-integrable systems in which all orbits are isolated. Then, 
an orbit must undergo a bifurcation when trMA = +2. 

The elements of the stability matrix Ma can be calculated from solutions of the linearized 
equation of motion in the transverse x direction, which we write in the Newtonian form 



Sx{t) = . (4) 

x=0,y=yA{t) 



This equation is identical to Hill's equation ||l^ in its standard form |11] 



6x{t) + [X + Q{t)]6x{t) = 0, (5) 

where Q{t) is a Ta (or Ta/2) periodic function whose constant Fourier component is zero. In 
general, the solutions of (^) are non-periodic. However, periodic solutions with period Ta (or Ta/2) 
and multiples thereof exist for specific discrete values of A. This happens exactly at bifurcations 
of the A orbit where trMA = +2. The periodic solutions 5x{t) found at these discrete values of A 
describe the x motion of the bifurcated orbits infinitely close to the bifurcation point. A special 
case of the Hill equation is the Lame equation which we discuss in the following section. 



3 The periodic Lame functions 



One of the standard forms of the Lame equation reads |12| 

A"{z) + \h- n{n + 1) sn^{z, k)] A{z) = , 



(6) 



where sn{z, k) is a Jacobi elliptic function with modulus k limited by < A; < 1. The real period of 
sn(z, k) in the variable z is 4K, where 



K = K{k)=F{-,k 



IT 



(7) 



is the complete elliptic integral of the first kind with modulus k. We follow throughout this paper 
the notation of Gradshteyn and Ryzhik [O for elliptic functions and integrals. We are interested 
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here only in real solutions for A(z) with real argument z. Hence h and n{n + 1) are assumed here 
to be arbitrary real constants. This means that n is either real, or complex with real part — ^. 
There is a vast literature on the periodic solutions of Eq. (^); for an exhaustive presentation of 
their definition and series expansions as well as the most relevant literature, we refer to Erdelyi et 
al |jl2|. (See also |Q for literature prior to 1932.) Ince ^, § has introduced a unique classification 
and nomenclature for the four types of periodic solutions, calling them Ec5^(2;) and Ec™(z), where 
n is the parameter appearing in (^) and m an integer giving the number of zeros in the interval 
< z < 2K. Following a slight redefinition by Erdelyi Q], the Ec(z) are even and the Es(2;) are 
odd functions of z — K, respectively. When m is an even integer, the Lame functions have the 
period 2K in the variable z, which is the same as the period of sn^(z,/c) appearing in (^); when 
m is odd, they have the period 4K. Solutions with period 2pK (p = 3, 4, . . .) can also be found; 
we shall discuss some solutions with period 8K further below. All these periodic solutions exist 
only for discrete eigenvalues of h, denoted by a™ and 6^ for the Ec™ and Es^, respectively; there 
exists exactly one solution of each of the above four types of Lame functions for each m > 0. The 
eigenvalues of h can, in principle, be found by solving the characteristic equation obtained from 
an infinite continued fraction which is, however, rather difficult to evaluate in the general case. 
In the context of our paper, they are determined by bifurcations of a linear periodic orbit and we 
obtain them therefore from a numerical calculation of its stability discriminant trMA. 

The Fourier expansions derived by Ince [Q, with the modification by Erdelyi |14], are given in 
terms of the variable 



c 



vr 



- - am(z. A;), 



where am(z, k) = arcsin[sn(2;, k)], and read as follows: 



Ecr(z) 
Ecr+i(z) 



Esr(z) 



Esr+i(z) 



-^o + ^^2r cos(2rC), 

^ r=l 

oo 

^yi2r+icos[(2r + l)C], 

r=0 

oo 

^52rSin(2rC), 



r=l 

oo 



^S2r+isin[(2r + l)C] 

r=0 



(period 2K) 
(period 4K) 
(period 2K) 
(period 4K) 



(8) 

(9) 
(10) 

(11) 
(12) 



with m = 0, 1, 2, . . . The expansion coefficients can be calculated by two-step recurrence relations; 
we give them here only for the 



[n{n + 1) - 2] k'^A2 = [2h - k^n{n + 1)]Aq , 



[n(n + 1) - (2r + 2)(2r + 1)] k^A 



2r+2 



2h - k^n{n + 1) - ^r'^{2 - k^) 



A 



2r 



- [n{n + 1) - (2r - 2)(2r - 1)] fc^A 



2r-2 , 



(13) 



(14) 



(with r = 1, 2, 3, . . . ) and refer to Erdelyi et al (|jT2|, ch 15.5.1) for the other recurrence relations 
which look very similar. 

Although the series (|9|) - (|l2|) are known to converge for A; < 1, they turned out to be 
semiconvergent in our numerical calculations for the cases with complex n, due to the fact that the 
characteristic values of h were only determined approximately. We have truncated the above series 
at the value Vmax where the corresponding coefficient has its smallest absolute value before starting 
to diverge. The cut-off values rmax were found to increase with the order m of the Lame function; 
their values are given in the Tables | and | in Sees. ^ and ^, respectively. 
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When n is an integer, the Fourier series terminate at finite values of r. The Lame functions then 
become finite polynomials in the Jacobi elliptic functions sn(z), dn(z), and cn(z), and are called 
the "Lame polynomials" in short. In Sec. ^ we will encounter a special case of the Lame equation 
in which h and n are not independent, but where h = ^n(n + 1) along with k"^ = ^. Then, to 
each integer n there exists only one value of m. Although the lowest few polynomials of this type 
and their eigenvalues of h are included in the tables given by Ince we give below their explicit 
expressions which take a particularly simple form. The basic four types of solutions correspond to 
the four rest classes modulo 4 of the integer n. With p = 0, 1, 2, 3, we obtain the following sums 
which are finite since the expansion coefficients become identically zero for r > p: 



EcJ(z) = 5]^4rcn4''(^), (period 2K) (15) 

P 

Ec|tl(^) = cn{z)J2CirCn^''{z), (period 4K) (16) 

r=0 

P 

EsJ+2(^) = dn(z)sn(z)^L'4rCn^''(^), (period 4K) (17) 

r=0 

P 

EsJ+3(z) = cn(z)dn(z)sn(2)^S4r.cn^''(z). (period 2K) (18) 

The simple one-step recurrence relations for the coefficients are (with r = 0, 1, 2, . . . , p — 1) 

(r+l)(4r + 3)^4r+4 = -[p(4p + 1) - r(4r + 1)] ^4r- , (19) 

2(r + l)(4r + 5)C74^+4 = - [(2p + l)(4p + 1) - (2r + l)(4r + 1)] C4,. , (20) 

2(r + l)(4r + 3)Z)4r+4 = -[(2p + l)(4p + 3) - (2r + l)(4r + 3)] 1)4,. , (21) 

(r + l)(4r + 5)^4^+4 = -[(p + l)(4p + 3) - (r + l)(4r + 3)] 54r • (22) 

The first sixteen Lame polynomials obtained from the above equations are given explicitly in Sec. 
P (see Table ^) , with the normalization Aq = Bq = Cq = Do = 1. 

For half-integer values of n, one obtains periodic Lame functions with period 8K that have 
algebraic forms in the Jacobi elliptic integrals and are called "algebraic Lame functions" ^ |l5| . For 
the special case with k"^ = ^ and h = ^n{n + 1) we will encounter them in Sec. |^ at period-doubling 
bifurcations of island-chain type. As shown in a beautiful paper by Ince , there exist two linearly 
independent periodic solutions for n = 2p + ^ with p = 0, 1,2, ... 

p p-i 



= \/dn(z) + cn{z) { sn^''{z) + cn{z) dn{z) J] Br sn'^''{z) \ , (23) 

r=0 r=0 
P P-1 



EC+i/2(^) = \/dn(z)-cn(z) { J] Ar sn^'\z) - cn{z) dn{z) J] B^ sn'^''{z) \ , (24) 



r=0 r=0 



where m is the number of zeros in the open interval (0, 2K); the coefficients A^ and Br are given 
by two coupled recurrence relations. For k"^ = ^, h = ^n{n + 1) = 2p{p + 1) + |, there is only one 
solution with m = p for each p, and the recurrence relations read 

2(2r + 2)2^^+1 + [Ap{p + 1) - 6r(2r + 1)] A - 4(p - r + l)(p + r) A-i 

= 2(2r + 2)Br+i - 3(2r + l)Br + 2rBr-i, (25) 
(2r + 2fBr+i + [2p{p + 1) - 3(r + l)(2r + 1)]S,, -2{p-r){p + r + l)Br-i 

= {2r + 2)Ar+i. (26) 
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These relations hold for r > provided that coefficients with negative indices r are taken to be 
zero. It is quite easy to see that Br = Ar+i = for r > p, which justifies the upper limits of the 
sums above. The coefficient Aq may be used for the overall normalization of both functions. 

The algebraic Lame functions ( p3| ) and (^) are even and odd functions of z, respectively, and 
related to each other by Ec^^^'^2('2 + 2K) = ^s^j^/^2^z), which amounts to a sign change in 
front of cn(z). Erdelyi [15| showed that the linear combinations Ec^^^^gC-^) + Es^^^'^2(-^) 

-^'^^+1/2 ('^) ~ ■^^2p+i/2('^) even and odd functions of z — K, respectively, and proposed that 
they be used instead of the functions (23,24) introduced by Ince. However, as we shall see in Sec. 

both pairs of independent functions are relevant in connection with period-doubling bifurcations. 
The Lame functions found for < p < 3 are explicitly given in Table ^ of Sec. ^ 



4 The Henon-Heiles potential 



We investigate here the role of the straight-line librating orbit A in the Henon-Heiles (HH) potential 
[l9|l. The Hamiltonian reads in scaled coordinates 



e = 6H = 6 



HH(x,y) = ^(x2 + y2)+^2y_ ly3^ (27) 



_2 

whereby the scaled energy is e = 1 at the saddle points. The Newton equations of motion are 

x + {l + 2y)x = 0, (28) 
y + y-y^ + x^ = 0. (29) 

These equations, and therefore the classical dynamics of the HH potential, depend only on the 
scaled energy e as a single parameter. In our numerical investigations we have solved Eqs. (28,29|) 



and determined the periodic orbits by a Newton- Raphson iteration using their stability matrix |lf:]. 

The basic periodic orbits with shortest periods found in the HH potential have been discussed 
mathematically by Churchill et al [^]; an exhaustive numerical search and classification has been 
performed by Davies et al We focus here on the straight-line A orbit which exists along the three 
symmetry axes of the HH potential, one of which coincides with the y axis. This orbit undergoes 
an infinite series of bifurcations which were studied in Ref . Q . They form a geometric progression 
on the scaled energy axes e, cumulating at the critical energy e = 1 where the period Ta becomes 
infinity and the orbit A becomes non-compact. All the orbits bifurcated from it exist, however, also 
at e > 1 and stay in a bounded region of the (x, y) space. Vieira and Ozorio de Almeida |jl^ have 
investigated some of these orbits at e > 1, both numerically and semi-analytically using Moser's 
converging normal forms near a harmonic saddle. In Fig. [l| we show the shapes of the orbits born at 
the isochronous bifurcations of orbit A, i.e., of those orbits having the same period Ta as orbit A at 
the bifurcation points. The subscripts of their names indicate their Maslov indices a needed in 



the context of the semiclassical periodic orbit theory |20| , 21, Although we make no use of the 
Maslov indices in the present paper, they are a convenient means of classification of the bifurcated 
orbits, as will become evident from the systematics below. 

All orbits shown in Fig. |l| are evaluated at the barrier energy e = 1. In the upper part of the 
figure, the x axis has been zoomed by a factor 0.163 from each panel to the next, in order to bring 
the shapes to the same scale. The orbits look practically identical in the lower 97% of their vertical 
range, but near the barrier (y = 1) they make one more oscillation in the x direction with each 
generation. In the lower part of the figure, we have zoomed also the y axis by the same factor from 
one panel to the next and plotted the top part of each orbit, starting from y = 1. In these blown-up 
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Figure 1: Orbits bifurcated from the A orbit in the HH potential, evaluated at energy e = 1. 
The subscripts give their Maslov indices a. Dashed lines: rotations Ro-, solid lines: librations 
(only one libration orbit is shown for each index; its partner is obtained by reflection at the vertical 
symmetry line containing orbit A). Top panels: successive scaling of x axis from left to right with 
the factor 0.163. Bottom panels: successive scaling of both axes with the same factor; along the y 
axis only the top part (starting from y = 1) is shown. 

scales, the tips of the orbits exhibit a perfect self-similarity which has been described by analytical 
scaling constants in Ref. [^. 

In Fig. |2| we show the stability discriminant trM for the orbit A (with its Maslov index a 
increasing by one unit at each bifurcation) and the orbits born at its isochronous bifurcations, 
plotted versus scaled energy e. In the lowest panel, we see the uppermost 3% of the energy scale 
available for the orbit A. The first bifurcation occurs at 65 = 0.969309, where A5 becomes unstable 
(with trMA > 2) and the stable orbit R5 is born. At eg = 0.986709, orbit Ag becomes stable 
again and a new unstable orbit Lg is born. In the middle panel, we have zoomed the uppermost 
3% of the previous energy scale. Here the behavior of A repeats itself, with the new orbits R7 and 
Lg born at the next two bifurcations. Zooming with the same factor to the top panel, we see the 
birth of Rg and Liq. This can be repeated ad infinitum: each new figure will be a replica of the 
previous one, with all the Maslov indices increased by two units and with trMA oscillating forever. 
This fractal behaviour is characteristic of the "Feigenbaum route to chaos" [|l|, |2|, ^ , although the 
present system is different from the Henon map in that the pitchfork bifurcations seen in Fig. ^ here 
are isochronous due to the reflection symmetry of the HH potential around the lines on which the 



bifurcating orbit A is situated (see Refs. |23, 24, |25[ for a discussion of the non-generic nature of the 
bifurcations in potentials with discrete symmetries). Also, the successive bifurcations happen here 
from one and the same orbit A, whereas in the standard Feigenbaum scenario one studies repeated 
period doubling bifurcations. 

Note that the functions trM(e) of the bifurcated orbits in Fig. |2| are approximately linear and 
intersect at e = 1 in two points, one for the librating orbits Lg- with trML(e = 1) = -1-8.183 (lying 
outside the figure), and one for the rotating orbits R^- with trMR(e = 1) = —4.183. We shall derive 
this linear behaviour in the limit e — > 1 from an asymptotic analytical evaluation of trMA in a 
forthcoming publication [p6|. 
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Figure 2: Trace of stability matrix M of orbit A and the orbits born at successive pitchfork 
bifurcations in the HH potential, plotted versus scaled energy e. From bottom to top: successively 
zoomed energy scale near e = 1. 

Presently we focus on the shapes of the orbits born at the bifurcations of A. Infinitely close to 
the bifurcation points, their motion in the transverse x direction is given by periodic solutions of 
the stability equation (|5[). Note that this equation is identical with the full equation of motion ( p^ 
in the x direction, which happens to be linear for the HH potential. The function i/Ait) describing 
the A orbit can easily be found analytically |Q and is given, with the initial condition yA(0) = yi, 
by 

yA{t) = yi + {y2-yi)sn^iz,k), (30) 
where z is the scaled time variable 



z = \l {ys - yi)/6 t = at, 



and yi {i = 1, 2, 3) are the roots of the cubic equation e = 6 Visi(,x = 0, y) 
are the turning points of the A orbit, whose period is 



(31) 

3y2 _ 2y3_ aj^^ 



Ta = 2^6/(y3 - yi) K = (2/a) K . (32) 

The modulus of the elliptic integral is given by 

fc^ = (2/2 - yi)/ (ys - yi) (33) 

and tends to unity for e — > 1 where 2/2 = 2/3- Rewriting Eq. ([2^ ) in terms of the scaled time variable 
z, it becomes identical with the Lame equation (P), with 

/i = 6(l + 2yi)/(y3-2/i), n(n + l) = -12 ^ n = -1/2 ± {i/2)V^. (34) 
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Table 1: Bifurcation energies e^, e* of the A orbit in the Henon-Heiles (HH) potential, names 
Oct of the bifurcated orbits, Lame functions of their motion Xa{t); cut-off order r^ax of the Fourier 
expansion, and period P of the Lame functions. The constant a is given in (31). Right part: 
isochronous bifurcations; left part: period-doubhng bifurcations. 



Note that h depends on the energy e. Hence, the discrete eigenvalues h = be directly 

related to the bifurcation energies eg- of the orbit A, and the corresponding Lame functions Ec™(z) 
and Es™(z) to the motion x{z) = x(at) of the new periodic orbits born at the bifurcations. In 
Table || we give the bifurcation energies obtained from the numerical computation of tr Ma, the 
names of the bifurcated orbits, and the corresponding Lame functions with their periods in the 
variable z = at given in (31). We also give the values r^ax at which the Fourier series (^) - ([l^ 
have been truncated. The right part of the table contains the lowest isochronous bifurcations seen 
in Fig. ^ and the bifurcated orbits shown in Fig. ||; their Lame functions all have the period 2K. 

The left part of Table |l] contains the lowest non-trivial period-doubling bifurcations (where 
trMA = —2), which are also of pitchfork type, and the names of the orbits born thereby. Their 
Lame functions have the period 4K. To avoid ambiguities, we denote their bifurcation energies by 
e*. The period-doubling bifurcations of new orbits with the Maslov indices 11, 12, 15, 16, etc., are 
trivial in the sense that they just involve the second iterates of orbit A and of the bifurcated orbits 
R5, Lg, R7, Lg, etc. The shapes of the first six non-trivial new orbits born at these bifurcations are 
shown in Fig. ^; they have similar scaling properties as those shown in Fig. |^. 
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Figure 3: The periodic orbits born at the six lowest period-doubling bifurcations of orbit A in the 
Henon-Heiles potential. The x axis has been scaled like in Fig. |l[ 
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Mathematically speaking, the periodic solutions for Xo-(i) in the form of the Lame functions exist 
only at the bifurcation energies However, the bifurcated orbits exist for all e > eo-- As long as 
the amplitude of their x motion remains small, it must be given by the Lame equation with the 
constants appearing in (|34[). But this equation has only periodic solutions when h has an eigenvalue 
corresponding to a bifurcation energy e^j. Therefore, the bifurcated orbits must keep their y motion 
"frozen" at y(t) = 2/^(0 with the parameters corresponding to Co-. Consequently, they also keep 
their periods at the bifurcation values. This has been confirmed numerically, as noticed already in 
iQl, to hold up to e = 1 and even beyond. Within the same small-amplitude limit of the x motion, 
the energy of the y motion is frozen at its value eo-, and the excess energy e — is consumed to 
rescale the amplitude of x{t). In other words, we can determine the normalization of the Lame 
function Xo-(i) of each bifurcated orbit by exploiting the energy conservation. This is most easily 
done at the time to = 7^/2 where y{€) has its maximum value, i.e., y(to) = V2-, y(io) = 0, and 
around which we know the symmetry of the Lame functions. For the even functions EcJ^ we have 
a^o-(*o) = 0) and Xo-(to) is, with (|27|), found to be 



(35) 



(36) 



Xo(to) = v/(e - e,)/3(l + 2y2) . 
For the odd functions Es™ we have Xo-(to) = 0, and their slopes at to ^'^^ given by 

i(T(io) = y^(e - e<x)/3 . 

In this way we can not only normalize the Lame functions near the bifurcation points, but also 
predict their evolution at higher energies. 

In Figs. § - show some of the periodic orbits obtained numerically from solving the equations of 
motion ( pstpsj ) at e = 1 by solid lines, and compare them to those predicted in the frozen-y-motion 
approximation, using y(t) = yAit) (given at their bifurcation energies or e*) and using for x{t) 
the Lame functions according to Table |l|, scaled as explained above. We see that in all cases, even 
for the lowest bifurcations, the new orbits keep their y motion acquired at their bifurcation energies, 
up to e = 1, indeed: the two curves y{t) and yA{t) can hardly be distinguished. As a consequence, 
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Figure 4: Orbit R5 in the HH potential at e = 1. Left panels: y{t) and x{t) versus time t (in units 
of 0.01); Right panel: orbit in the (x, y) plane. Solid lines: numerical results obtained by solving 
Eqs. (|2^,29). Dashed lines: y{t) given by yA{t) in (^0|) at the bifurcation energy 65, and x{t) given 
by theXame function according to Table scaled as described in the text. 
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Figure 5: Same as Fig. ^ for the orbit Fg. Note that this orbit has twice the period of orbit A 
at eg and is given by the Lame function Ec^(at) with period 4K and m = 3 zeros in the interval 
< ai < 2K. 



we can expect the functions x{t) to be well described by the appropriate Lame functions. This is, 
indeed, the case if the latter are correctly scaled. As we see, the normalization predicted by (35,36D 
is the better, the closer the bifurcation energy comes to the saddle-point energy e = 1. A rigorous 
justification of the frozen-y-motion approximation will be presented elsewhere j26[ |. 

We point out that all orbits born at the isochronous pitchfork bifurcations in the HH system are 
given by Lame functions with period 2K, since the orbit A is given by yA{t) ( |30|) and hence has the 
same period as the function appearing in the Lame equation. The Lame functions with period 4K 
must therefore correspond to orbits born at period-doubling bifurcations (see Table |l] and Fig. ^). 



0.008 



^ 0.0 



-0.008 




1.0 



0.5 



0.0 



-0.5 




250 500 750 1000 1250 
tx 100 



-0.007 0.0 0.007 
X 



Figure 6: Same as Fig. |^ for the orbit Lg. 
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Figure 7: Same as Fig. ^ for the orbit R13. Note that x{t) is harmonic over a long time while y is 
close to the saddle {y — 1); the period of this harmonic oscillation is uj±_ = 27rAy/3, as shown in Ref. 
[Q] and also obtained from an asymptotic expansion of the Lame functions [p6| . 



Note also that according to bifurcation theory |23, 25, 27 1, two degenerate periodic orbits 
should be born at each isochronous pitchfork bifurcation. The librating orbits Lo- come, indeed, 
in pairs that are symmetric to the y axis, whereas the rotating orbits Ro- can be run through in 
two opposite directions. Each of these pairs of orbits are, however, described by one and the same 
Lame function for x{t) which is invariant under the corresponding symmetry operation. This is 
in agreement with a theorem, proved by Ince [^, that there cannot exist two linearly independent 
periodic solutions of the Lame equation to the same characteristic value of h. An exception of this 
theorem is given by transcendental and algebraic Lame functions with integer or half-integer values 
of n (see Sec. ^ for an example). 

We finally note that in the figures § - |^, the time scales are not always normalized such that 
t = Q corresponds to y(0) = yi, as assumed in Eq. (|30|), but obtained rather randomly due to the 
way in which the periodic orbits where searched and found numerically. However, if we shift the 
time origin to t' = according to (pO[), all figures illustrate how the Lame functions Ec^(2;) are 
even and the Es™(z) odd, according to their definition [|^, |T^, around t' = K/a where y{t') has its 
maximum. This demonstrates that the association of Lame functions to bifurcated orbits allows 
one to understand (or predict) their symmetries. 



5 The quartic Henon-Heiles potential 

We next investigate the quartic Henon-Heiles (H4) potential [p^, p9[ with the scaled Hamiltonian 



AH 



i;{pl+pl) + ynA{x,y) 



^H4(x, y) = \ (x2 + y^)-\ [x^ + y') + \ xW, (37) 



which is similar to the HH potential but has four saddles at the scaled energy e = 1; it has reflection 
symmetry at both coordinate axes and both diagonals. It contains straight-line librating orbits 
along all four symmetry lines; two of them, which we call again orbits A, oscillate between the two 
saddles lying on the coordinate axes. To be specific, we choose again the A orbit along the y axis. 
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It has the same behaviour as the A orbit in the HH potential, but it approaches a saddle at both 
ends. Its motion is, for 2/a(0) = 0, given by 

VAit) = yisn{at,k) , a = y2/V2, k = yi/y2, (38) 

and its period is 

Ta = 4^/2K/2/2 = 4K/a. (39) 
Hereby ±yi and ±y2 are the solutions of e = 4VH4(a; = 0,y) = 2?/^ — y^, i.e., 

yi = ^1- VT^, y2 = \/lWT^, (40) 

and ibyi are the turning points of the orbit. 

The linearized equation of motion in the x direction, which decides about the stability of the 
orbit A, is for the H4 potential 

x(t) + [l + 3y2(t)]x(t) =0, (41) 

neglecting here explicitly a term of order x^. Inserting the solution for yA{t) in (|3^ ) and transforming 
to the scaled time variable z = at leads again to the Lame equation (^ with 

h = 2/yl, n(n+l) = -6 ^ n = -1/2 ± (i/2)\/23 . (42) 

Compared to the HH potential, we have now a new situation which is a consequence of the higher 
symmetry of the H4 potential: the periodic function sv?{z,k) appearing in the stability equation 
has half the period, namely 2K, of that of the orbit A itself ([39|). Therefore, all its periodic solutions 
with period 2K, corresponding to Lame functions with an even number m of zeros, also have the 
period Ta/2 at the bifurcations. The solutions involving the Lame functions with odd m share their 
periods 4K = oTa with that of the A orbit. 

The systematics of the isochronous bifurcations of the A orbit for increasing bifurcation energies 
Co- IS given m Table H The new orbits appear with m = 2, 3, 4, . . . , with alternatingly odd and 
even Lame functions. Like in the HH case, the EcJ^ correspond to librations L^- and the EsJ^ to 
rotations R^-. They appear alternatingly as 2K and 4K periodic functions. The orbits given in the 
left part of the table are born stable and remain stable up to e > 1, whereas those in the right part 
are unstable at all energies. 
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Table 2: The same as in Table |^, but only for the isochronous bifurcations of the A orbit in 
the quartic Henon-Heiles (H4) potential. Note the alternating appearance of 2K and 4K periodic 
Lame functions with increasing bifurcation energy eg-. The constant a is given in Eq. (|38|). Orbits 
appearing on the left side are stable up to e > 1, those on the right side are unstable at all energies. 



12 



10 




1 1 1 


— 1 1 1 1 






8 












6 








He/ 


1-10 I-'g / 


S 4 
2 








M 




As 






Ay/ \ \ Ag 





-2 








/ ^\ 


y 


1 . 1 . 1 . 1 . 1 .1.1.1 



I I I I I I I I I I I I I I I 

0.85 0.875 0.9 0.925 0.95 0.975 1.0 1.025 

e 



Figure 8: Stability discriminant of the shortest orbits in the H4 potential, plotted versus scaled 
energy e (cf. Tabled for the bifurcation energies). Lj bifurcates at e ~ 1.027 and becomes Lg. 

The non-trivial period-doubling bifurcations in this potential are of island-chain type, and the 
orbits born thereby are given by Lame functions of period 8K. We will not investigate them here, 
but refer to the analogous situation in the quartic oscillator potential discussed in Sec. |6|. 

In Fig. ^ we show the stability discriminant tr M versus energy e for the orbit A and the orbits 
born at its lowest isochronous bifurcations. Different from the HH potential, here the functions 
trM(e) of the bifurcated orbits are, to a good approximation, quadratic in e. This can be derived 
analytically ||2^. It is also striking that, different from Fig. |2|, trM of all the orbits shown is always 
larger than or equal to —2. This behaviour, together with the systematics seen in Table ^, can be 
explained by the following arguments. 

Bearing in mind that the stability matrix of the second iterate of a periodic orbit O is just 
Mq, where Mq is that of the primitive orbit, one easily finds that its discriminant is 

tr Mo2 = tr = (tr Mq)^ - 2 , (43) 

which can never be less than —2. Hence, we can mathematically understand trM of the orbit A in 
Fig. ^ to be that of a second iterate. Its primitive is half of the orbit A, having the same period 
as that of the function sn^ in the Lame equation for its stability, and having a discriminant tr M 
which oscillates around zero, exceeding the values +2 and —2 on both sides symmetrically, like 
trMA in the HH potential (Fig. |2|). Their second iterates, which correspond to the full bifurcated 
orbits, therefore have a discriminant trM which is quadratic in e. These features, together with the 
systematics in Table are all a consequence of the C^^ symmetry of the H4 potential, including the 
reflection symmetry at the x axis that divides the A orbit (and all the bifurcated orbits discussed 
here) into to equal halves. The quadratic behaviour of tr M of the bifurcated orbits is also consistent 
with the fact that their next period-doubling bifurcations (where tr M = —2) are symmetry breaking 
(see Ref. [H for details). 

Like in the HH case, the values of trM of the bifurcated orbits intersect in two points at e = 1, 
one with trM(e = 1) = —1.711 for the orbits born stable, and one with trM(e = 1) = +9.991 for 
the orbits born unstable. 

We thus obtain the result that Lame functions with both periods 2K and 4K describe the orbits 
born at isochronous bifurcations of the A orbit. Two examples are shown in Figs. |^ and |lO|, where 
the X motion of the orbits L14 and Rie is given by the functions Ec^(2;) and Es^(2:), respectively. 
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Figure 9: Same as Fig. ^ for the orbit L14 in the H4 potential (cf. Table ^). Note that x{t) has 
half the period of y{t). 

Whereas the latter shares its period with orbit A, the former has half its period. Like before, these 
orbits are evaluated at the critical energy e = 1. The normalization of the Lame functions has been 
chosen as for the HH potential, using the frozen-y-motion approximation and energy conservation, 
leading here to 

= ^J{e-e^)/2{l + 3yi) = ^{e-e„)/2. (44) 
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Figure 10: Same as Fig. H for the orbit Rie in the H4 potential (cf. Table 
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6 The homogeneous quartic oscillator 



We now turn to the quartic oscillator (Q4) potential 

^Q4(x,y) = Vx^ + /)+ixV (45) 



which has been the object of several classical and semiclassical studies [30|. Since it is homogeneous 
in the coordinates, the Hamiltonian can be rescaled together with coordinates and time such that 
its classical mechanics is independent of energy. Consequently, the system parameter is here not 
the energy but the parameter e. The potential (|45| ) has the same symmetry as the H4 potential ( |37|) 
and, correspondingly, possesses periodic straight-line orbits along both axes. The motion of the A 
orbit along the y axis is given by 

yA{t) = yocn{yot,k), yQ = {AEfl\ ^ = 1/2 (46) 

with the period Ta = 4K/?/o- Its turning points are ±?/o- Note that this solution does not depend 
on the value of e. The stability of the orbit A, however, does depend on e. The linearized equation 
of motion for the transverse x motion yields, after transformation to the coordinate z = yot, the 
Hill equation 

x"{z) + €[l-sii^{z,k)]x{z) =0. (47) 
This is a special case of the Lame equation with 

h = e = ^n{n + l), (48) 

where we have used k'^ = 1/2. The nice feature is that we here know analytically the eigenvalues 
h = hn of the Lame equation, namely those given in Eq. (Esl). This agrees with the analytical result 



for the stability discriminant of the A orbit (|4q), which has been derived long ago by Yoshida |31]: 

tr Ma = 4 cos 



+ 2 . (49) 



It is easy to see that the bifurcation condition trMA = +2 leads exactly to the values ( |48| ) of the 
parameter e. 

The periodic solutions of the Lame equation (^) at the bifurcation values of /i = e are the Lame 
polynomials discussed already in Sec. |3[ Their explicit forms for n = 0, 1, . . . , 15 are given in Table 
^, using the short notation cn = cn{z, k), etc., and a normalization such that their leading coefficient 
is unity. We also give in Table S the names of the new-born orbits, using the same nomenclature as 



for the H4 potential in the previous section. Their shapes are shown in Fig. 11. They have exactly 
the same topologies as the orbits of the H4 potential and are again described by Lame functions of 
pairwise alternating periods 2K and 4K, as seen also in Table ^. Each of these orbits has a discrete 
degeneracy of 2, which is due to the time reversal symmetry for the rotations and to the reflection 
symmetries about the coordinate axes for the librations. 

A special comment is due concerning the cases n = 1 and n = 2 which correspond to e = 1 



and 3, respectively. For these values of the parameter e, the Q4 potential is integrable 31] 
and no bifurcations occur for the shortest orbits. Nevertheless, the Lame equation ( p7|) possesses 
mathematically the solutions Ec} and Es2, respectively. The orbits B3 and C4 given in Table ^ and 



Fig. [11| have topologically the shapes given by these Lame polynomials, but it should be emphasized 
that they are not generated through bifurcations but are generic orbits existing at all values of e. 
Under the symmetry operation e — > (3 — e)/(l -|- e), which corresponds to a rotation about 45 



degrees and simultaneous stretching of the potential [M, the orbits of type A are mapped onto the 
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Table 3: Numbers n and parameter values e„ for the isochronous bifurcations of the A orbit in the 
homogeneous quartic oscillator (Q4) potential; names Oq- of bifurcated orbits; Lame polynomials 
for their x motion, and their periods P. The Jacobi elliptic functions sn, cn, and dn are given in 
short notation as sn = sn(z,fc), etc., with z = yot, where yo and the modulus k are given in (46). 
(See text for a special comment on orbits B3 and C4.) 



orbits of type B and vice versa, and the orbits of type C are mapped onto themselves. This is seen 
easily in the Lame polynomial describing the B orbit, EcJ(2;) = cn{z,k), which is proportional to 
the function (^) describing the A orbit. 



L3 [B3] [C4] R; L7 Rg Rp 




Figure 11: Shapes of orbits born at the isochronous bifurcations of orbit A in the Q4 potential. 
Vertical direction: yA{t) given by (^6[); horizontal direction: x{t) given by the Lame polynomials 
listed in Table ^ (with scaled normalizations to fit the shapes into the frames). 
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The scaling properties of these orbits and their evolution away from the bifurcation values e^^ 
are more difficult to analyze than in the HH and H4 potentials and will be discussed elsewhere p^ . 

We now discuss period-doubling bifurcations of the A orbit in the Q4 potential. There is a series 
of trivial period doublings which just involve the second iterates of the orbit A and the orbits born 
at its isochronous bifurcations and listed in Table ^. The non-trivial period doublings occur when 
trMA = — 2, which leads with (H^) to the critical values 



e = 2p(p + l)+3/8, P = 0,l,2, 



(50) 



This is exactly one of the conditions Q for the existence of period-8K solutions of (^7|), namely that 
obtained by inserting n = {4p + l)/2 into Eq. (^8|). The solutions are the algebraic Lame functions, 
given (up to p = 3) in Table Q below. These bifurcations are of the island-chain type (see, e.g., Ref. 
pTf ) : the quantity tr M of the second iterate of orbit A - let us call it orbit A^ - touches the value 
+2, but the orbit A^ remains stable on either side. At the bifurcation, two doubly-degenerate orbits 



are born, one stable and one unstable. The situation is illustrated in Fig. 12 around the bifurcation 
at e = 4 + 3/8 {p = 1). The unstable new orbit is here called Fio, and the stable new orbit is called 
Pg. There shapes, together with those born at the other period doublings listed in Table ^, are 
shown in Fig. |l^ below. Their degenerate symmetry partners are called F^ for the librating orbits 
(obtained by reflecting the orbits at the x axis) and for the rotating orbits (obtained by time 
reversal of the Po- orbits). 

According to the theory of Ince Q, the algebraic Lame functions of period 8K are one exceptional 
case where two independent periodic solutions can coexist for the same critical value of h. These are 
the functions Fc2p^y2 -^^2p+i/2 defined in Eqs. ( p3| , p^ ). As we see from Table |^, they correspond 
to the unstable orbits of type F„ and F^. In contrast to the degenerate pairs of bifurcated orbits 
in the HH and H4 potentials, which are represented by one and the same periodic Lame function, 
the pairs Fg- and F^ are here given by two linearly independent functions. With this, however, 
the number of independent solutions of the second-order differential equation (^) is exhausted. 
Therefore the other pair of stable orbits of type Pq- and P^ born at the period doublings cannot be 
given by any new independent solutions. Indeed, we see from Table § that the orbits P^ and P^ are 





Figure 12: Stability discriminant trM of period-2 orbits in the Q4 potential around a period- 
doubling bifurcation of island-chain type at e = 4.375. A^ is the second iterate of the primitive 
orbit A5 (shown by a dashed line); its Maslov index increases from 9 to 11 at the bifurcation. Fio 
and Pg are the genuine period-2 orbits bifurcating from it; each has a discrete degeneracy of two 
(see text and Table §); Pg bifurcates at e ~ 5.3 and becomes Ps. 
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Table 4: Numbers n and parameter values e„ of the four lowest non-trivial period-doubling bifur- 
cations of the A orbit in the homogeneous quartic oscillator. Given are also the names Oo- of the 
four topologically different orbits born simultaneously at these bifurcations, and the (linear combi- 
nations of) algebraic Lame functions with period 8K that describe their x motion. (The same short 
notation for the Jacobi elliptic functions sn, dn and cn is used as in Table |.) 



given by the two independent linear combinations Ec2^_^(^2 ^ -^^2^+1/2 '^^i^h were constructed by 
Erdelyi to have the same symmetry properties as those of the 2K and 4K periodic Lame 
functions, as discussed in Sec. ^. 

We thus have found the interesting result - which was new to us - that the stable and unsta- 
ble pairs of orbits born at a period-doubling bifurcation of island-chain type are mutually linear 
combinations of each other. It follows from our above arguments that this result must hold for all 
Hamiltonians with the double reflection symmetry C2V 
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Figure 13: Same as Fig. O, but for orbits born at period-doubling bifurcations of orbit A in the 
Q4 potential. The motion x{t) is given by the algebraic Lame functions listed in Table ^. 
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Figure 14: The four orbits born at the period-doubling bifurcation at e = 12.375. Left: their 
shapes in the {x,y) plane; right: algebraic Lame functions describing their motion x{t). 

In Fig. lU we illustrate the situation for the four orbits born at the bifurcation at e = 12 + 3/8 
{p = 2). In the left panels, their shapes x{y) are shown; in the right panels we plot the algebraic 
Lame functions (and their linear combinations) which describe the motion x(t) of the respective 
orbits. Looking merely at the shapes of these orbits, the pairwise linear dependence of their x 
motion is not obvious at all. 



7 Summary and conclusions 

We have investigated cascades of isochronous pitchfork bifurcations of straight-line librational orbits 
in two-dimensional potentials. The linearized equation of the x motion transverse to these orbits, 
determining their stability, can be written in the form of the Lame equation. Its eigenvalues cor- 
respond to the bifurcation values of the system parameter, given also by the condition trM= +2 
for the stability discriminant of the straigh-line orbit, and its eigenfunctions describe the x motion 
of the new orbits born at the bifurcations. These eigenfunctions are the periodic Lame functions 
of period 2K or 4K, where K is the complete elliptic integral determining the period of the parent 
orbit at the corresponding bifurcation. In potentials with C2v symmetry, the solutions occur alter- 
natingly as Lame functions of period 2K and 4K, respectively. When this symmetry is absent, the 
4K periodic solutions describe the orbits born at period-doubling pitchfork bifurcations. 

We have shown numerically that the periodic Lame functions describe very accurately the shapes 
of the bifurcated orbits obtained from a numerical integration of the equations of motion, as long as 
the amplitude of their x motion remains small, i.e., as long as one is not too far from the bifurcation 
point. Exploiting the energy conservation in the Henon-Heiles type potentials HH and 114 and the 
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known symmetries of the Lame functions, we can predict the propagation of the new orbits up 
to the critical saddle-point energy where they have all become unstable and the system is highly 
chaotic. We thus have found an analytical description of an infinite series of unstable periodic orbits 
in chaotic systems. 

In the homogeneous quartic oscillator (Q4) potential, the series expansions of the periodic Lame 
functions terminate and they become finite polynomials. In this potential we have also analyzed 
solutions of period 8K which occur at period-doubling bifurcations of the straight-line orbits of 
island-chain type. The two pairs of orbits born thereby are represented by two independent sets 
of orthogonal periodic solutions of the Lame equation, which here are identified with the so-called 
algebraic Lame functions that can again be given in a closed form. 

Similar cascades of pitchfork bifurcations have also been discussed in connection with the dia- 



magnetic Kepler problem represented by hydrogen atoms in strong magnetic fields |24, |3^, 34]. 
Expressing the Hamiltonian in (scaled) semiparabolic coordinates {u,v), the effective potential for 
orbits with angular momentum = (where the z axis is the direction of the external mag- 
netic field) becomes similar to the H4 and Q4 potentials discussed here (although it contains only 
quadratic and sixth-order terms in the coordinates). Since physically the {u,v) coordinates are pos- 
itive definite, the periodic orbits in the diamagnetic Kepler problem correspond to the half-orbits 
of the 114 and the Q4 potentials. There exist straight-line librating orbits, corresponding to oscilla- 
tions of the electron along the symmetry axis, which bifurcate infinitely many times as the energy 
of the electron approaches the ionization threshold. The stability of these linear orbits is given by 
the Mathieu equation, which is analogous to the Lame equation (|6|) but with the function sn'^(z, k) 
replaced by cos(2z). Its periodic solutions are the periodic Mathieu functions se^ and cem which 
have properties completely analogous to those of the periodic Lame functions, and were actually 
studied in detail by Ince [35| prior to his investigations of the Lame functions. The topology of the 
Mathieu functions and of the bifurcated orbits described by them is exactly the same as for the 
Q4 and 114 potentials described here. In particular, the so-called "balloon" orbits B„ and "snake" 
orbits S„ with n = 1,2, ... correspond exactly to the alternating sequence of halves of the 
orbits R5, Rg, . . . and L7, Ln, . . . shown in Fig. |l^. (The other orbits, born unstable at the bifur- 
cations, were not considered in |24| since they do not pass through the centre.) We believe that our 
analysis, applied in terms of the Mathieu functions, may be useful for further investigations of the 
bifurcations occurring in the diamagnetic Kepler problem. 

The knowledge of the analytical properties of the bifurcated orbits will be useful in the appli- 
cation of the periodic orbit theory to the potentials studied here. First steps in this direction have 
been quite successful [28, but the orbits bifurcated from the A orbit were not considered. 

Their incorporation into the semiclassical trace formula is the object of further work in progress. 
We expect the bifurcated orbits, in particular, to play an important role in the semiclassical analysis 
of resonances above the barriers in Henon-Heiles type or similar potentials. 



Acknowledgments 

We are grateful to S Fedotkin and A Magner for very stimulating discussions and critical comments. 
Valuable comments by M Sieber and H Then are highly appreciated. We also acknowledge financial 
support by the Deutsche Forschungsgemeinschaft. 



20 



References 



[1] Feigenbaum M J 1978 J. Stat. Phys. 19 25 
see also Feigenbaum M J 1983 Physica 7 D 16 

Bountis T C 1981 Physica 3 D 577 

Greene J M, McKay R S, Vivaldi F and Feigenbaum M J 1981 Physica 3 D 577 

Brack M 2001 Festschrift in honor of the 75th birthday of Martin Gutzwiller ed A Inomata et 
al Foundations of Physics 31 209 [LANL preprint |nlin.CD/000603^ 

Lame G 1839 Journal de Mathematiques Pures et Appliquees (Liouville) 4 126 
Lame G 1839 ibid. 4 351 

see, e.g., Strutt M J O 1932 Ergebnisse der Mathematik und ihrer Grenzgehiete Vol. 1 no. 3 
(Springer- Verlag, Berlin) 

Ince E L 1940 Proc. Royal Soc. Edinburgh 60 47 
Ince E L 1940 Proc. Royal Soc. Edinburgh 60 83 
Henon M and Heiles C 1964 Astr. J. 69 73 
Hill G W 1886 Acta Math. 8 1 

Magnus W and Winkler S 1966 Hill's Equation (Interscience Publ., New York) 

Erdelyi A et al 1955 Higher Transcendental Functions Vol HI (McGraw-Hill, New York) ch 15. 
Beware of the misprint in equation (13) and elsewhere in chapter 15.5.1 of this reference; the 
correct definition of the variable C is that given in Eq. (P) of our present paper 

Gradshteyn I S and Ryzhik I M 1994 Table of Integrals, Series, and Products (Academic Press, 
New York, 5th edition) ch 8.1 

Erdelyi A 1941 Phil. Mag. (7) 31 123 

Erdelyi A 1941 Phil. Mag. (7) 32 348 

Tanaka K and Brack M, to be published 

Churchill R C, Pecelli G and Rod D L 1979 Stochastic Behavior in Classical and Quantum 
Hamiltonian Systems ed G Casati and J Ford (Springer- Verlag, New York) p 76 

Davies K T R, Huston T E and Baranger M 1992 Chaos 2 215 

Vieira W M and Ozorio de Almeida A M 1996 Physica 90 D 9 

Gutzwiller M C 1971 J. Math. Phys. 12 343 

Gutzwiller M C 1990 Chaos in classical and quantum mechanics (Springer, New York) 

Brack M and Bhaduri R K 1997 Semiclassical Physics, Frontiers in Physics Vol. 96 (Addison- 
Wesley, Reading, USA) 

de Aguiar MAM, Malta C P, Baranger M, and Davies K T R 1987 Ann. Phys. (NY) 180 167 

Mao J-M and Delos J B 1992 Phys. Rev. A 45 1746 

Then H L 1999 Diploma thesis (Universitat Ulm) 
Then H L and Sieber M, to be published 

Fedotkin S N, Magner A G, Mehta M and Brack M, to be published 

Sieber M 1996 J. Phys. A 29 4715 
Schomerus H and Sieber M 1997 J. Phys. A 30 4537 
Sieber M and Schomerus H 1998 J. Phys. A 31 165 

Brack M, Creagh S C and Law J 1998 Phys. Rev. A 57 788 

Brack M, Meier P and Tanaka K 1999 J. Phys. A 32 331 



21 



[30] Eckardt B 1988 Phys. Rep. 163 205 

Bohigas O, Tomsovic S and Ullmo U 1993 Phys. Rep. 233 45 
Eriksson A B and Dahlqvist P 1993 Phys. Rev. E 47 1002 

Lakhshminarayan A, Santhanam M S, and Sheorey V B 1996 Phys. Rev. Lett. 76 396 

[31] Yoshida H 1984 Celest. Mech. 32 73. Strictly speaking, Yoshida derived trM for the half-orbit 
to be trMA/2 = 2\/2cos[7rVl + 8e/4], from which Eq. (|4|) follows using Eq. (|3|) 

[32] Mehta M and Brack M, work in progress 

[33] Main J, Wiebusch G, Holle A and Welge K H 1986 Phys. Rev. Lett. 57 2789 
Main J, Wiebusch G, Holle A and Welge K H 1987 Z. Phys. D 6 295 

[34] Wintgen D and Friedrich F 1987 Phys. Rev. A 36 131 
Friedrich F and Wintgen D 1989 Phys. Rep. 183 39 

Hasegawa H, Robnik M and Wunner G 1989 Prog. Th. Phys. Suppl. 98 198 
[35] Ince E L 1932 Proa. Royal Soc. Edinburgh 52 355 

[36] Brack M, Bhaduri R K, Law J and Murthy M V N 1993 Phys. Rev. Lett. 70 568 

Brack M, Bhaduri R K, Law J, Maier Ch and Murthy M V N 1995 Chaos 5 317,707(E) 



22 



